/**
 * @author Gao Jian
 * @email gaoj@njupt.edu.cn
 * @create date 2021-06-03 10:23:40
 * @desc
 * 关于tiff掩膜数据的处理，包括连接破碎的小块成大块，删除块内空洞，较小的块
 */

#define cimg_use_tiff

#include <CImg.h>
#include <cuda_runtime.h>
#include <gdal/gdal.h>
#include <iostream>

using namespace std;
using namespace cimg_library;

// void make_histogram(unsigned *hist_i, unsigned *hist_o, unsigned char *data,
//                     unsigned char *mark, unsigned nh, unsigned nd);
void make_vind(unsigned char *datao, unsigned char *datav, unsigned char *dataw,
               unsigned char *datab, unsigned char *datar, float ratio,
               unsigned sv, unsigned sw, unsigned sb, unsigned sz,
               unsigned char verbose);
void PrintDeviceInfo();

int main(int argc, char **argv) {
  const char *fr =
      cimg_option("-fr", (char *)(NULL), "Reference TiFF File Name");
  const char *fv = cimg_option("-fv", (char *)(NULL), "NDVI TiFF File Name");
  const char *fw = cimg_option("-fw", (char *)(NULL), "NDWI TiFF File Name");
  const char *fb = cimg_option("-fb", (char *)(NULL), "NDBI TiFF File Name");
  const char *fo = cimg_option("-fo", (char *)(NULL), "Output TiFF File Name");

  const char disp = cimg_option("-dis", 0, "Visualize data");
  const char *pmode =
      cimg_option("-p", "D9J9F0E9C0", "Procedure mode, 'DJFEC'");

  const float ratio = cimg_option("-rt", 1., "ratio.");

  const unsigned verbose = cimg_option("-v", 0, "verbose.");
  const unsigned device_id = cimg_option("-d", 0, "device ID.");

  int deviceCount;
  cudaGetDeviceCount(&deviceCount);

  if (verbose) {
    PrintDeviceInfo();
    printf("使用GPU设备: %d.\n", device_id);
  }
  if (device_id >= deviceCount) {
    if (verbose)
      cout << "错误的设备ID " << device_id << ", 应小于 " << deviceCount
           << "!!!" << endl;
    return -1;
  }
  cudaSetDevice(device_id);

  if (verbose)
    printf("========== IS. ==========\n");

  if (!fr || !fv || !fw || !fb)
    return -1;

  // =========================================================================
  GDALAllRegister();
  GDALDatasetH hDatasetR = GDALOpen(fr, GA_ReadOnly);
  GDALDatasetH hDatasetV = GDALOpen(fv, GA_ReadOnly);
  GDALDatasetH hDatasetW = GDALOpen(fw, GA_ReadOnly);
  GDALDatasetH hDatasetB = GDALOpen(fb, GA_ReadOnly);
  GDALDatasetH hDatasetO = GDALOpen(fo, GA_Update);
  if (!hDatasetR || !hDatasetV || !hDatasetW || !hDatasetB || !hDatasetO)
    return -1;

  //   double adfGeoTransform[6];

  //   GDALGetGeoTransform(hDataset, adfGeoTransform);
  //   unsigned resolution = resl;
  //   double arc = 6378245. * 3.14159265358 * 2 / 360;
  //   if (resl == 0) {
  //     if (strstr(GDALGetProjectionRef(hDataset), "degree"))
  //       resolution =
  //           unsigned((adfGeoTransform[1] - adfGeoTransform[5]) / 2 * arc +
  //           0.5);
  //     else
  //       resolution =
  //           unsigned((adfGeoTransform[1] - adfGeoTransform[5]) / 2 + 0.5);
  //   }

  // int nw, nh;
  GDALRasterBandH hBandR = GDALGetRasterBand(hDatasetR, 1);
  GDALRasterBandH hBandV = GDALGetRasterBand(hDatasetV, 1);
  GDALRasterBandH hBandW = GDALGetRasterBand(hDatasetW, 1);
  GDALRasterBandH hBandB = GDALGetRasterBand(hDatasetB, 1);
  GDALRasterBandH hBandO = GDALGetRasterBand(hDatasetO, 1);
  int nw = GDALGetRasterBandXSize(hBandR);
  int nh = GDALGetRasterBandYSize(hBandR);
  // GDALGetBlockSize(hBand, &nw, &nh);
  // printf(
  //     "Block=%dx%d Type=%s, ColorInterp=%s\n", nw, nh,
  //     GDALGetDataTypeName(GDALGetRasterDataType(hBand)),
  //     GDALGetColorInterpretationName(GDALGetRasterColorInterpretation(hBand)));
  CImg<unsigned char> datar(nw, nh);
  CImg<unsigned char> datao(nw, nh);
  CImg<unsigned char> datav(nw, nh);
  CImg<unsigned char> dataw(nw, nh);
  CImg<unsigned char> datab(nw, nh);
  CPLErr error;
  error = GDALRasterIO(hBandR, GF_Read, 0, 0, nw, nh, datar.data(), nw, nh,
                       GDT_Byte, 0, 0);
  error = GDALRasterIO(hBandO, GF_Read, 0, 0, nw, nh, datao.data(), nw, nh,
                       GDT_Byte, 0, 0);
  error = GDALRasterIO(hBandV, GF_Read, 0, 0, nw, nh, datav.data(), nw, nh,
                       GDT_Byte, 0, 0);
  error = GDALRasterIO(hBandW, GF_Read, 0, 0, nw, nh, dataw.data(), nw, nh,
                       GDT_Byte, 0, 0);
  error = GDALRasterIO(hBandB, GF_Read, 0, 0, nw, nh, datab.data(), nw, nh,
                       GDT_Byte, 0, 0);

  datav -= datav.min();
  dataw -= dataw.min();
  datab -= datab.min();
  datav /= 2;
  dataw /= 2;
  datab /= 2;

  int ntv = datav.max();
  int ntw = dataw.max();
  int ntb = datab.max();

  // =========================================================================
  if (verbose)
    printf("========== Dimension [%d] [%d] [%d]. ==========\n", ntv + 1,
           ntw + 1, ntb + 1);

  // CImg<unsigned int> vv_i(ntv + 1, 1, 1, 1, 0);
  // CImg<unsigned int> vw_i(ntw + 1, 1, 1, 1, 0);
  // CImg<unsigned int> vb_i(ntb + 1, 1, 1, 1, 0);

  // CImg<unsigned int> vv_o(ntv + 1, 1, 1, 1, 0);
  // CImg<unsigned int> vw_o(ntw + 1, 1, 1, 1, 0);
  // CImg<unsigned int> vb_o(ntb + 1, 1, 1, 1, 0);

  // make_histogram(vv_i.data(), vv_o.data(), datav.data(), datar.data(), ntv +
  // 1,
  //                nw * nh);
  // make_histogram(vw_i.data(), vw_o.data(), dataw.data(), datar.data(), ntw +
  // 1,
  //                nw * nh);
  // make_histogram(vb_i.data(), vb_o.data(), datab.data(), datar.data(), ntb +
  // 1,
  //                nw * nh);

  // CImg<float> vind(ntv, ntw, ntb, 1, 0.);

  // make_vind(datao.data(), datav.data(), dataw.data(), datab.data(),
  //           datar.data(), ratio, ntv + 1, ntw + 1, ntb + 1, nw * nh,
  //           verbose);

  int nw_t = 3;
  int nh_t = 3;

  CImg<unsigned char> datar_t(nw_t, nh_t, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0);
  CImg<unsigned char> datao_t(nw_t, nh_t);
  CImg<unsigned char> datav_t(nw_t, nh_t, 1, 1, 2, 1, 0, 2, 1, 0, 1, 0, 0);
  CImg<unsigned char> dataw_t(nw_t, nh_t, 1, 1, 2, 1, 0, 2, 1, 0, 1, 0, 0);
  CImg<unsigned char> datab_t(nw_t, nh_t, 1, 1, 2, 1, 0, 2, 1, 0, 1, 0, 0);

  int ntv_t = datav_t.max();
  int ntw_t = dataw_t.max();
  int ntb_t = datab_t.max();

  make_vind(datao_t.data(), datav_t.data(), dataw_t.data(), datab_t.data(),
            datar_t.data(), ratio, ntv_t + 1, ntw_t + 1, ntb_t + 1, nw_t *
            nh_t, verbose);

  // =========================================================================
  error = GDALRasterIO(hBandO, GF_Write, 0, 0, nw, nh, datao.data(), nw, nh,
                       GDT_Byte, 0, 0);

  GDALClose(hDatasetR);
  GDALClose(hDatasetV);
  GDALClose(hDatasetW);
  GDALClose(hDatasetB);
  GDALClose(hDatasetO);

  // =========================================================================

  if (verbose)
    printf("========================= <<<< Finished. >>>> "
           "=========================\n");
  // if (fo)
  //   dataout.save(fo);

  if (disp)
    (datar, datao).display("REFERENCE | OUTPUT");

  return 0;
}
